Zimbabwe Soybean Trials

Loading data

#install.packages("sf")
#install.packages("stars")
library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
data <- read.csv("data/Zimbabwe.csv",header = T,sep = ";")
# colnames(data)
# head(data)
pheno<- data[, c(4, 7, 16,21,25,27,38,41,46,63,64,65,66,68)] 
table(pheno$LOCATION_NAME)
## 
##      Banket_SCL     Bindura_SCL Chisumbanje_CBI   Goromonzi_SCL       Gwebi_CBI 
##             288             288             270             558             198 
##      Harare_CBI      Harare_KKR      Harare_SCL      Kadoma_CBI      Mazowe_TCK 
##             468              90             288             180              36 
##   Mthampden_TCK      Mutare_SCL      Shamva_CBI  Stapleford_SCL 
##              90             378             360             468
colnames(pheno)
##  [1] "CROP_SEASON"     "TRIAL_NAME"      "LOCATION_NAME"   "REP_NO"         
##  [5] "Crop_Variety"    "Control"         "PL_Lodging"      "R8_PL_Height"   
##  [9] "GRAIN_YLD_13PCT" "PLANT_DATE"      "HARV_DATE"       "SITE_LAT"       
## [13] "SITE_LONG"       "IRRIGATION"
pheno$geno <- gsub(" ", "_",pheno$Crop_Variety)
pheno$geno <- factor(pheno$geno)
# unique(pheno$IRRIGATION)
levels(pheno$geno) <- paste("G", 1:length(levels(pheno$geno)), sep = "_")
# table(pheno$Crop_Variety)
# table(pheno$geno)

pheno[which(pheno$HARV_DATE==""),"LOCATION_NAME"]
##   [1] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##   [5] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##   [9] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [13] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [17] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [21] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [25] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [29] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [33] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [37] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [41] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [45] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [49] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [53] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [57] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [61] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [65] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [69] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [73] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [77] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [81] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [85] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [89] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [93] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
##  [97] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
## [101] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
## [105] "Banket_SCL"     "Banket_SCL"     "Banket_SCL"     "Banket_SCL"    
## [109] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [113] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [117] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [121] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [125] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [129] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [133] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [137] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [141] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [145] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [149] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [153] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [157] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [161] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [165] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [169] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [173] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [177] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [181] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [185] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [189] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [193] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [197] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [201] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [205] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [209] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [213] "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"    "Bindura_SCL"   
## [217] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [221] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [225] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [229] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [233] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [237] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [241] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [245] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [249] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [253] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [257] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [261] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [265] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [269] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [273] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [277] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [281] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [285] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [289] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [293] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [297] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [301] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [305] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [309] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [313] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [317] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [321] "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL"  "Goromonzi_SCL" 
## [325] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [329] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [333] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [337] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [341] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [345] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [349] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [353] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [357] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [361] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [365] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [369] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [373] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [377] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [381] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [385] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [389] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [393] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [397] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [401] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [405] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [409] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [413] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [417] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [421] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [425] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [429] "Harare_SCL"     "Harare_SCL"     "Harare_SCL"     "Harare_SCL"    
## [433] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [437] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [441] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [445] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [449] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [453] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [457] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [461] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [465] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [469] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [473] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [477] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [481] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [485] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [489] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [493] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [497] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [501] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [505] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [509] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [513] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [517] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [521] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [525] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [529] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [533] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [537] "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"     "Mutare_SCL"    
## [541] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [545] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [549] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [553] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [557] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [561] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [565] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [569] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [573] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [577] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [581] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [585] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [589] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [593] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [597] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [601] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [605] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [609] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [613] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [617] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [621] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [625] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [629] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [633] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [637] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [641] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [645] "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL" "Stapleford_SCL"
## [649] "Harare_CBI"
pheno$DTH <- as.Date(as.character(pheno$HARV_DATE), format="%m/%d/%Y")-
                  as.Date(as.character(pheno$PLANT_DATE), format="%m/%d/%Y")-5
# dim(pheno)
# colnames(pheno)
# head(pheno)
# str(pheno)
pheno$GRAIN_YLD_13PCT <-  gsub(",", ".", pheno$GRAIN_YLD_13PCT)
trial <-  gsub("_.*", "", pheno$LOCATION_NAME)
pheno$trial <- paste0(trial,"_", pheno$CROP_SEASON)
unique(pheno$TRIAL_NAME) # 41 trials
##  [1] "S18_2_ZW_SCL5"  "S18_2_ZW_SCL11" "S18_2_ZW_CBI4"  "S18_2_ZW_SCL7" 
##  [5] "S18_2_ZW_CBI1"  "S18_2_ZW_SCL6"  "S18_2_ZW_CBI2"  "S18_2_ZW_SCL10"
##  [9] "S18_2_ZW_CBI9"  "S18_2_ZW_SCL8"  "S19_1_ZW_CBI2"  "S19_1_ZW_CBI4" 
## [13] "S19_2_ZW_SCL5"  "S19_2_ZW_CBI4"  "S19_2_ZW_SCL1"  "S19_2_ZW_CBI1" 
## [17] "S19_2_ZW_SCL2"  "S19_2_ZW_SCL6"  "S19_2_ZW_CBI5"  "S19_2_ZW_SCL3" 
## [21] "S20_2_ZW_SCL10" "S20_2_ZW_SCL11" "S20_2_ZW_SCL7"  "S20_2_ZW_CBI2" 
## [25] "S20_2_ZW_CBI1"  "S20_2_ZW_SCL8"  "S20_2_ZW_SCL12" "S20_2_ZW_SCL9" 
## [29] "S21_2_ZW_SCL3"  "S21_2_ZW_SCL2"  "S21_2_ZW_CBI2"  "S21_2_ZW_CBI1" 
## [33] "S21_2_ZW_KKR1"  "S21_2_ZW_CBI3"  "S21_2_ZW_SCL1"  "S22_2_ZW_SCL3" 
## [37] "S22_2_ZW_CBI1"  "S22_2_ZW_TCK1"  "S22_2_ZW_CBI3"  "S22_2_ZW_SCL4" 
## [41] "S23_2_ZW_SCL1"  "S23_2_ZW_TCK1"  "S23_2_ZW_SCL2"
unique(pheno$LOCATION_NAME) # 14 Locations
##  [1] "Banket_SCL"      "Bindura_SCL"     "Chisumbanje_CBI" "Goromonzi_SCL"  
##  [5] "Harare_CBI"      "Harare_SCL"      "Kadoma_CBI"      "Mutare_SCL"     
##  [9] "Shamva_CBI"      "Stapleford_SCL"  "Gwebi_CBI"       "Harare_KKR"     
## [13] "Mazowe_TCK"      "Mthampden_TCK"
unique(pheno$CROP_SEASON) # 6 years
## [1] 2018 2019 2020 2021 2022 2023
length(unique(pheno$Crop_Variety)) #98 genotypes
## [1] 98
head(pheno)
##   CROP_SEASON    TRIAL_NAME LOCATION_NAME REP_NO  Crop_Variety Control
## 1        2018 S18_2_ZW_SCL5    Banket_SCL      1     S1079/6/7      No
## 2        2018 S18_2_ZW_SCL5    Banket_SCL      1  TGx 1987-62F      No
## 3        2018 S18_2_ZW_SCL5    Banket_SCL      1  TGx 1991-22F      No
## 4        2018 S18_2_ZW_SCL5    Banket_SCL      1       Lukanga      No
## 5        2018 S18_2_ZW_SCL5    Banket_SCL      1    S1146/5/25      No
## 6        2018 S18_2_ZW_SCL5    Banket_SCL      1 TGx 2014-16FM      No
##   PL_Lodging R8_PL_Height GRAIN_YLD_13PCT PLANT_DATE HARV_DATE SITE_LAT
## 1          1          105         4680.55 12/04/2018 4/26/2019 -176,425
## 2          5           95         2236.97 12/04/2018 4/26/2019 -176,425
## 3          5           90         1866.21 12/04/2018 4/26/2019 -176,425
## 4          1          100         3723.68 12/04/2018 4/26/2019 -176,425
## 5          1           90         3185.52 12/04/2018 4/26/2019 -176,425
## 6          5          100         2851.49 12/04/2018 4/26/2019 -176,425
##    SITE_LONG IRRIGATION geno      DTH       trial
## 1 30,399,997    Rainfed G_25 138 days Banket_2018
## 2 30,399,997    Rainfed G_56 138 days Banket_2018
## 3 30,399,997    Rainfed G_57 138 days Banket_2018
## 4 30,399,997    Rainfed  G_8 138 days Banket_2018
## 5 30,399,997    Rainfed G_27 138 days Banket_2018
## 6 30,399,997    Rainfed G_69 138 days Banket_2018
pheno <- pheno |> rename("Plant Height"=8, "Grain Yield"=9, "Plant lodging"=7, "rep"=4)

traits <- c(colnames(pheno)[c(7,8,9)])
traits
## [1] "Plant lodging" "Plant Height"  "Grain Yield"
pheno[, traits] <- lapply(pheno[, traits], as.numeric)
## Warning in lapply(pheno[, traits], as.numeric): NAs introduced by coercion
factors <- c(colnames(pheno)[c(1,2,3,4,5,6,14,15,16,17)])
factors
##  [1] "CROP_SEASON"   "TRIAL_NAME"    "LOCATION_NAME" "rep"          
##  [5] "Crop_Variety"  "Control"       "IRRIGATION"    "geno"         
##  [9] "DTH"           "trial"
traits <- c(colnames(pheno)[c(7,8,9)])
traits
## [1] "Plant lodging" "Plant Height"  "Grain Yield"
pheno[, factors] <- lapply(pheno[, factors], as.factor)
head(pheno)
##   CROP_SEASON    TRIAL_NAME LOCATION_NAME rep  Crop_Variety Control
## 1        2018 S18_2_ZW_SCL5    Banket_SCL   1     S1079/6/7      No
## 2        2018 S18_2_ZW_SCL5    Banket_SCL   1  TGx 1987-62F      No
## 3        2018 S18_2_ZW_SCL5    Banket_SCL   1  TGx 1991-22F      No
## 4        2018 S18_2_ZW_SCL5    Banket_SCL   1       Lukanga      No
## 5        2018 S18_2_ZW_SCL5    Banket_SCL   1    S1146/5/25      No
## 6        2018 S18_2_ZW_SCL5    Banket_SCL   1 TGx 2014-16FM      No
##   Plant lodging Plant Height Grain Yield PLANT_DATE HARV_DATE SITE_LAT
## 1             1          105     4680.55 12/04/2018 4/26/2019 -176,425
## 2             5           95     2236.97 12/04/2018 4/26/2019 -176,425
## 3             5           90     1866.21 12/04/2018 4/26/2019 -176,425
## 4             1          100     3723.68 12/04/2018 4/26/2019 -176,425
## 5             1           90     3185.52 12/04/2018 4/26/2019 -176,425
## 6             5          100     2851.49 12/04/2018 4/26/2019 -176,425
##    SITE_LONG IRRIGATION geno DTH       trial
## 1 30,399,997    Rainfed G_25 138 Banket_2018
## 2 30,399,997    Rainfed G_56 138 Banket_2018
## 3 30,399,997    Rainfed G_57 138 Banket_2018
## 4 30,399,997    Rainfed  G_8 138 Banket_2018
## 5 30,399,997    Rainfed G_27 138 Banket_2018
## 6 30,399,997    Rainfed G_69 138 Banket_2018
#Irrigation system
# table(pheno$trial,pheno$IRRIGATION, useNA = "ifany")
# unique(pheno$IRRIGATION)
pheno[which(pheno$trial=="Banket_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Banket_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Bindura_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Chisumbanje_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Goromonzi_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Goromonzi_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Gwebi_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Mutare_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial==" Mutare_2020"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Stapleford_2019"),"IRRIGATION"] <- "Rainfed"
pheno[which(pheno$trial=="Stapleford_2020"),"IRRIGATION"] <- "Rainfed"

# table(pheno$IRRIGATION)
# 
# 
# levels(pheno$IRRIGATION)

Data visualization

pheno_long <- reshape2::melt(pheno, measure.vars = traits, variable.name = "trait") #transformando o pheno em formato longo
# head(pheno_long)

ggplot(pheno_long, aes(x = trait, y = value, color = trait)) + geom_boxplot(outlier.shape = NA)+
    geom_point(alpha = 0.2, position = "jitter", size = 1) + facet_grid(trait~LOCATION_NAME,
    scales = "free_y") + theme_bw() + theme(axis.text.x = element_text(angle=90))
## Warning: Removed 191 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 191 rows containing missing values or values outside the scale range
## (`geom_point()`).

library(ggh4x)
library(ggpmisc)
## Loading required package: ggpp
## Registered S3 methods overwritten by 'ggpp':
##   method                  from   
##   heightDetails.titleGrob ggplot2
##   widthDetails.titleGrob  ggplot2
## 
## Attaching package: 'ggpp'
## The following object is masked from 'package:ggh4x':
## 
##     stat_centroid
## The following object is masked from 'package:ggplot2':
## 
##     annotate
data_vis=ggplot(data = pheno_long, aes(x=trait, y= value, color=trait)) +
  geom_boxplot() +
  geom_point( position = "jitter", alpha = 0.5) +  facet_nested(trait~LOCATION_NAME+CROP_SEASON, scales = "free_y") +
  theme_bw();data_vis
## Warning: Removed 191 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Removed 191 rows containing missing values or values outside the scale range
## (`geom_point()`).

locations <-  pheno[,c("TRIAL_NAME","LOCATION_NAME","SITE_LAT","SITE_LONG")]
locations_unq <- locations |> distinct(SITE_LAT,.keep_all = TRUE)

locations_unq <- locations_unq[-21,]
smart_clean_coords <- function(x, decimal_places = 5) {
  # Replace commas with dots (fix decimals)
  x <- gsub(",", ".", x)
  
  # Remove thousand separators (dots between digits)
  x <- gsub("(?<=\\d)\\.(?=\\d)", "", x, perl = TRUE)
  
  # Convert to numeric (handles scientific notation like 0.00E+00)
  num <- as.numeric(x)
  
  # Function to scale down numbers that are too large
  scale_if_needed <- function(val) {
    if (is.na(val)) return(NA)
    
    # Latitude should be between -90 and 90
    if (abs(val) > 90) {
      # Keep dividing by 10 until in valid range
      while (abs(val) > 90) {
        val <- val / 10
      }
      return(val)
    } else {
      return(val)  # Already valid, return as-is
    }
  }
  
  # Apply scaling only where needed
  num <- sapply(num, scale_if_needed)
  
  # Round to specified decimal places
  round(num, decimal_places)
}
locations_unq$SITE_LAT <- smart_clean_coords(locations_unq$SITE_LAT)
locations_unq$SITE_LONG <- smart_clean_coords(locations_unq$SITE_LONG)

write.csv(locations_unq,file = "data/loc_unq.csv",row.names = F)

Mapping the trials

#Com ggplot2

# elevação

library(raster)
library(RColorBrewer)
library(ggplot2)
library(sf)

africa.sf<-st_read("data/afr_g2014_2013_0.shp")
# uf.sf<- brasil.sf[brasil.sf$UF %in% c('MT','GO','RS','PE','BA','ES','DF',
#             'RJ','AL','SE','MG','MS','SP','PR','SC'), ]

library(stringi)
africa.sf$ADM0_NAME <- stri_replace_all_fixed(
  africa.sf$ADM0_NAME,
   "C\xf4te d'Ivoire",  # Literal special character
  "Cote d'Ivoire"
)
# elevação
#Selecionando os estados

africa<- shapefile('data/afr_g2014_2013_0.shp')
africa
#brasil[LINHA , COLUNA]


uf<- africa[africa@data$ADM0_NAME%in% c("Zimbabwe"),] #Mesma indexacao de data.frame



# # uf
# 
 plot(uf, col='orange')
# plot(uf, col=c("red","blue"))
# text(uf,'UF', cex=1, halo=T)

#Cortando um raster para o estado
agro<- raster('data/agro.tif') # Importando o Raster de elevação
hwsd <- raster('data/hwsd.tif')

plot(agro)
plot(hwsd)
#for( i in 1:5){
hwsd_africa<- crop(africa,uf)
plot(hwsd_africa)

hwsd_africa<-mask(africa,uf)
plot(hwsd_africa)
palete.spec<- brewer.pal("Spectral",n=11)
#plot(elev_uf, col=rev(palete.spec))

writeRaster(hwsd,filename= "data/hwsd_map", overwrite=TRUE)
writeRaster(agro,filename= "data/agro_map", overwrite=TRUE)
WGS84<-CRS("+proj=longlat +datum=WGS84 +no_defs ") #
#pontos
pontos<- read.csv("data/loc_unq.csv", header=TRUE, sep=',') # se estiver usando outro OS
str(pontos)
# pontos[which(pontos$latitude < -23),]
# pontos[which(pontos$Instituicao=="CATI/Avare"),"regiao"]="SE"
# pontos[which(pontos$Instituicao=="DETEC/Itai"),"regiao"]="SE"
# pontos[which(pontos$cidade=="Campos Novos"),"regiao"]="SU"

# Apply to your data (5 decimal rounding)
# pontos$SITE_LAT <- smart_clean_coords(pontos$SITE_LAT)
# pontos$SITE_LONG <- smart_clean_coords(pontos$SITE_LONG)

# Check results
head(pontos)
tail(pontos)

coordinates(pontos)<- ~ SITE_LONG + SITE_LAT
str(pontos)

#plot(pontos,pch=20,cex=0.5, add=T)
crs(pontos)<-WGS84


# fv <- unique(pontos$regiao) 
# (m <- match(pontos$regiao, fv)) 


### Plotando o mapa do Brasil 
#

ggplot() + 
  geom_sf(data = africa.sf)
  #scale_fill_discrete(type = c("blue",'yellow','red','magenta','green'))


### Plotando o mapa do Brasil  e da regiao cortada
#

ggplot() + 
  geom_sf(data = africa.sf, fill="white")


pontos.sf
pontos

### Aproximando para a regi?o cortada
#


pontos@coords.nrs
(bbox_uf<-st_bbox(uf)) #Pegando os limites da regi?o cortada
ggplot() + 
  geom_sf(data = africa.sf, fill="white") +
  #Código para definir os limites do mapa cood_sf()
    coord_sf(xlim =c(bbox_uf["xmin"],bbox_uf["xmax"]) ,ylim=c(bbox_uf["ymin"],bbox_uf["ymax"]))+
   theme_light()
  # 
  


### Adicionando os rótulos de estado
#

ggplot() + 
  geom_sf(data = africa.sf, fill="grey") +
  theme_light()
 

### Plotando pontos
#

(pontos.sf<-data.frame(pontos)) #retornando nossos pontos para tabela

# pontos.sf <-pontos.sf[-21,] 
pontos.sf
# ggplot() + 
#   geom_sf(data = brasil.sf, fill="grey") +
#   geom_sf(data = uf.sf,fill="white")+
#   coord_sf(xlim =c(bbox_uf["xmin"],bbox_uf["xmax"]) ,ylim=c(bbox_uf["ymin"],bbox_uf["ymax"])) +
#   geom_sf_label(data = brasil.sf, aes(label = UF),fill=NA,label.size = NA)+
#   #Código para plotagem de pontos (geom_point)
#   geom_point(data = pontos.sf, mapping = aes(x = longitude, y = latitude, color=regiao),size =1) + 
#   scale_color_manual(breaks = c("NE", "CO", "SE","SU"),values=rainbow(4))+
#   theme_light()+ guides(color = guide_legend(title = "Region"))+labs(y="Latitude",x="Longitude")


problems <- africa.sf$ADM0_NAME[!stringi::stri_enc_isascii(africa.sf$ADM0_NAME)]
if(length(problems) > 0) print(problems)


### com rasters
#
# install.packages("tmap")
library(stars)

elev_uf.stars<- read_stars("data/hwsd_map.grd",RasterIO = list(nBufXSize = 600, nBufYSize = 600))
# elev_uf.stars<- read_stars("data/agro_map.grd",RasterIO = list(nBufXSize = 600, nBufYSize = 600))
names(elev_uf.stars)<-"Elevation"
pontos

st_bbox(africa)
(bbox_uf<-st_bbox(uf))
(pontos.sf<-data.frame(pontos)) 



mapa=ggplot() + 
  geom_sf(data = africa.sf, fill="white") +
  #Código para plotagem do rater geom_stars()+
  coord_sf(xlim =c(bbox_uf["xmin"],bbox_uf["xmax"]) ,ylim=c(bbox_uf["ymin"],bbox_uf["ymax"]),crs = "WGS84") +theme_light()+
  geom_point(data = pontos.sf, mapping = aes(x = SITE_LONG, y = SITE_LAT),size =1,show.legend = T)+
    geom_sf_label(data = africa.sf, aes(label = ADM0_NAME),fill=NA,label.size = NA)+theme(legend.position = "none",axis.text.x =element_text(angle = 90))+labs(x="",y="");mapa


africamap=ggplot() + 
  geom_sf(data = africa.sf, fill="white")+
    scale_fill_distiller(palette = "RdYlGn",na.value = NA,name = "Altitude")+
  geom_stars(data = elev_uf.stars,alpha =0.8)+
  coord_sf(xlim =c(-25.35875,63.50265) ,ylim=c(-46.9813,  37.56095),crs = "WGS84")+labs(y="Latitude",x="Longitude")+
  theme_light()+theme(legend.position = "none")+
  geom_rect(data=africa.sf,
    xmin = 25.23703,
    ymin =-22.41774,
    xmax = 33.05630,
    ymax =-15.60883 ,
    fill = NA, 
    colour = "black",
    linewidth = 0.8
  );africamap


library(cowplot)

map=africamap |> ggdraw()+
  draw_plot(
    {
      mapa + 
        coord_sf(
          xlim = c(25.23703, 33.05630),
          ylim = c(-22.41774, -15.60883),
          expand = FALSE) +
        theme(legend.position = "none")
      },
    # The distance along a (0,1) x-axis to draw the left edge of the plot
    x = 0.13, 
    # The distance along a (0,1) y-axis to draw the bottom edge of the plot
    y = 0.05,
    # The width and height of the plot expressed as proportion of the entire ggdraw object
    width = 0.40, 
    height = 0.40);
map









ggsave(plot=map,device = "pdf",filename ="Plots/soybean_zimbabwe.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=map,device = "png",filename ="Plots/soybean_zimbabwe.png",dpi = "retina",height = 6,width = 8)

# Single-environment analyses

These analyses had a purpose: to test the significance of the genetic variance in each environment. We used the following model:

\[ \mathbf{y} = \mathbf{X_1b} + \mathbf{X_2che} + \mathbf{Z_1g} + \mathbf{\epsilon} \]

where \(\mathbf{y}\) is the vector of phenotypic observations, \(\mathbf{b}\) is the vector of fixed effects of replication, \(\mathbf{che}\) is the vector of fixed effects of checks, \(\mathbf{g}\) is the vector of random effects of genotypes and \(\mathbf{\epsilon}\) the vector of errors. \(\mathbf{X_1}\), \(\mathbf{X_2}\), \(\mathbf{Z_1}\) $ are incidence matrix of \(\mathbf{b}\), \(\mathbf{g}\) effects respectively.

At each environment, we filtered out trials that do not have significance for genotype effect. The Likelihood Ratio Test (\(LRT\)) was computed as follows:

\[ LRT= −2 \times (Log𝐿 - Log L_𝑅) \]

where \(L\) is the maximum point of residual likelihood function of the complete model and \(L_R\) is the same for the reduced model, that is, without the effect to be tested. The LRT value was compared with a tabulated value based on the chi-square table, with one degree of freedom and 0.95 probability.

We calculated the broad-sense heritability (\(h^2\)):

\[ h^2 = \sigma^2g / \sigma^2g + \sigma^2e \]

where \(\sigma^2g\) is the is the additive genetic variance,\(\sigma^2g\) is non additive genetic variance and \(\sigma_e^2\) is the error variance.

STMT

source("https://raw.githubusercontent.com/saulo-chaves/May_b_useful/main/update_asreml.R")
#source("GEBV/codes/outliers.R")

library(asreml)

data.list = split(pheno, f = pheno$trial)

traits
stmt = list()
j="Banket_2018" 
#j="2021NAC_ST01"
i="Plant lodging"
unique(pheno$trial)


for (j in names(data.list)) {
  
   x = data.list[[j]]
  cat("====> Environment:", j, fill = TRUE)
  


 #   A.env = full2sparse(K = amat[which(rownames(amat) %in%
 #                                              c(x$geno)),
 #                                     which(colnames(amat) %in%
 #                                            c(x$geno))])
 # attr(A.env, "INVERSE") = FALSE
 
#  H.env.van = full2sparse(K = H.inv.van[which(rownames(H.inv.van) %in%
#                                               c(x$geno)),
#                                      which(colnames(H.inv.van) %in%
#                                             c(x$geno))])
#  
# attr(H.env.van, "INVERSE") = TRUE


  # x = x |> mutate(nad=geno)  |> mutate_at(vars(env, rep, geno, row, col, check,nad), as.factor)
  
  summar.mat = matrix(nrow = length(traits), ncol = 4, 
                      dimnames = list(traits, 
                                      c("geno","error","H2", "CV")))

  vc.env = list()
  Blue.env = list()
  
  for (i in traits) {
    
    cat("====> Trait:", i, fill = TRUE)
    
    if(all(is.na(x[,i]))) next
    
    model.r = asreml(fixed = x[,i] ~ rep, 
                   random = ~ geno,maxit = 50,
                   data = x, na.action = na.method(x = "include", y = "include"))
    

    
   # colocar nugget and spline
    model.f = asreml(fixed = x[,i] ~ rep+geno, 
                    maxit = 50,
                   residual = ~units,
                   data = x, na.action = na.method(x = "include", y = "include"))
 
    if(!model.r$converge) next
    
    if(any(na.exclude(model.r$vparameters.pc > 1))) model.r = update.asreml(model.r)
    # varcomp = summary(modelhvan)$varcomp
    varcomp = summary(model.r)$varcomp
      varcomp.df = data.frame(
      effect = c('genotypic', "residual"),
      component = c(round(
        varcomp[grep('geno', rownames(varcomp)),1],4),
        round(varcomp[grep('!R', rownames(varcomp)),1],4))
    )
     # if(!modelad$converge) next
     # 
     # if(any(na.exclude(modelad$vparameters.pc > 1))) modelad = up.mod(modelad)
     # 
    if(!model.f$converge) next
    
    if(any(na.exclude(model.f$vparameters.pc > 1))) model = up.mod(model.f)
      
    lrtest = as.data.frame(lrt(asreml(fixed = x[,i] ~ rep, 
                        residual = ~units, maxit = 50,
                        data = x),model.r))
    lrtest=up.mod(lrtest)
      
    lrtest$effect = c("genotypic")
    
    varcomp= full_join(varcomp.df, lrtest[,-1], by = "effect")
    varcomp$signi = ifelse(varcomp$`Pr(Chisq)` <= 0.06, "*", "ns")
    varcomp$signi2 = ifelse(varcomp$`Pr(Chisq)` <= 0.06, 
                            paste(round(varcomp$component, 4), "*"), 
                            paste(round(varcomp$component, 4), "ns"))
    rownames(varcomp)=varcomp$effect
    
    pred = predict(model.r, classify = "geno", sed = TRUE)
    pred$sed= pred$sed[-which(pred$pvals$status=="Aliased"),-which(pred$pvals$status=="Aliased")]
    MVdelta = mean((pred$sed^2)[upper.tri(pred$sed^2,diag = F)])
    H2 = 1-(MVdelta/(2*varcomp[varcomp$effect == "genotypic","component"]))
    h2 = vpredict(model.r, h2~(V1)/(V1+V2))

    #h2ad = vpredict(modelad, h2~(V3)/(V3+V4))
 
    CV = sqrt(varcomp[grep("r", rownames(varcomp)),"component"])/
      mean(x[,i], na.rm = TRUE)*100
    
    summar.mat[i,] = c(varcomp$signi2[1], round(varcomp$component[2],4), 
                       round(h2$Estimate,4), round(CV,4))
    
    
    pred=predict(model.f, classify = "geno", vcov = T)

    pred$vcov= pred$vcov[-which(pred$pvals$status=="Aliased"),-which(pred$pvals$status=="Aliased")]
    pred$pvals= pred$pvals[-which(pred$pvals$status=="Aliased"),]

    pred$pvals$w=diag(solve(pred$vcov))
    pred=pred$pvals
    
  vc.env[[i]] = varcomp.df

  Blue.env[[i]] = pred
  }
  stmt[[j]] = list(Blues = Blue.env,summary = summar.mat,vc=vc.env)
  
  rm(summar.mat,GEBV.env, varcomp, lrtest, model, 
     h2, CV, x, H.inv.env)
}

saveRDS(stmt, file = "Saves/stmt.rds")
#rm(data.list)

#hist(data[which(data$env=="UGN2021ABI_ST02"),"rytha"])
#table(data[which(data$env=="UGN2021ABI_ST02"),"rytha"],useNA = "ifany")

#table(data[which(data$env=="UGN2020NAS_ST04"),"vir2"], useNA = "ifany")

Results

stmt=readRDS("Saves/stmt.rds")

stvc=lapply(stmt, function(x)as.data.frame(cbind(x$summary)))

#stvc.df=lapply(stmti, function(x)as.data.frame(cbind(x$vc)))



#stmt[[1]][["Blues"]][[1]][,-c(3,4)]


stvc=do.call(rbind,stvc)
stvc$trait=rownames(stvc)
stvc$env=sub("\\..*","",stvc$trait)
stvc$trait = sub(".*\\.","",stvc$trait)
stvc$geno_sig=sub("^\\S+\\s+","",stvc$geno)
stvc$geno=sub("\\s.*","",stvc$geno)
stvc$geno= as.numeric(stvc$geno)
stvc$H2=as.numeric(stvc$H2)
stvc$error=as.numeric(stvc$error)
stvc$CV=as.numeric(stvc$CV)
library(tidyverse)
head(stvc)
##                                  geno       error     H2      CV         trait
## Banket_2018.Plant lodging      0.7079      1.1989 0.3713 46.7037 Plant lodging
## Banket_2018.Plant Height     484.7230    123.1807 0.7974 10.5836  Plant Height
## Banket_2018.Grain Yield   165411.8384 324929.2135 0.3373 22.4875   Grain Yield
## Banket_2019.Plant lodging          NA          NA     NA      NA Plant lodging
## Banket_2019.Plant Height           NA          NA     NA      NA  Plant Height
## Banket_2019.Grain Yield   184876.1880 295453.6310 0.3849 17.2950   Grain Yield
##                                   env geno_sig
## Banket_2018.Plant lodging Banket_2018        *
## Banket_2018.Plant Height  Banket_2018        *
## Banket_2018.Grain Yield   Banket_2018        *
## Banket_2019.Plant lodging Banket_2019     <NA>
## Banket_2019.Plant Height  Banket_2019     <NA>
## Banket_2019.Grain Yield   Banket_2019        *
stvcl= stvc |> select(geno, error, trait, env, geno_sig) |> pivot_longer(geno:error) 
  

unique(stvc$env)
##  [1] "Banket_2018"      "Banket_2019"      "Banket_2020"      "Bindura_2018"    
##  [5] "Bindura_2020"     "Bindura_2021"     "Chisumbanje_2018" "Chisumbanje_2019"
##  [9] "Goromonzi_2018"   "Goromonzi_2019"   "Goromonzi_2020"   "Goromonzi_2021"  
## [13] "Goromonzi_2022"   "Goromonzi_2023"   "Gwebi_2020"       "Gwebi_2021"      
## [17] "Harare_2018"      "Harare_2019"      "Harare_2020"      "Harare_2021"     
## [21] "Harare_2022"      "Kadoma_2018"      "Kadoma_2021"      "Mazowe_2022"     
## [25] "Mthampden_2023"   "Mutare_2018"      "Mutare_2019"      "Mutare_2020"     
## [29] "Mutare_2021"      "Shamva_2018"      "Shamva_2019"      "Shamva_2022"     
## [33] "Stapleford_2018"  "Stapleford_2019"  "Stapleford_2020"  "Stapleford_2022" 
## [37] "Stapleford_2023"
# cols.label=c("A20NAC", "A20NAS", "B20NAS", "A21ABI" ,"B21ABI","A21NAC", "B21NAC", "A21NAS", "B21NAS")
# names(cols.label)=c(unique(stvc$env)) 

stvcl|>  ggplot(aes(x=name,y=value))+ geom_col(data= subset(stvcl, !grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(fill="*"), color="black", fill="#a8ddb5") + geom_col(data= subset(stvcl,  grepl("error",stvcl$name)),color="black",fill="lightblue") + facet_grid(trait~env,scales = "free_y") + geom_text(data=subset(stvcl, !grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(label = geno_sig), size=6 ) + theme(legend.position = "top", )  + geom_col(data= subset(stvcl, grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(fill="ns"), fill="gray", color="black")  + theme(legend.position = "top", ) + labs(x="Effects", y="Genotypic variance component") +   geom_text(data=subset(stvcl, grepl("ns", stvcl$geno_sig) & grepl("geno",stvcl$name)),aes(label = "ns", vjust = -1.0), size=3 ) 
## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Removed 17 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_text()`).

stvc_bef <- stvc |>  select(H2, CV, trait, env, geno_sig) |> pivot_longer(H2:CV) |>  
  ggplot(aes(x=env,y=value))+geom_col(aes(fill = trait), position = position_dodge())+ facet_grid(name~trait, scales = "free_y") + theme(legend.position = "none", axis.text.x = element_text(angle = 90)) + labs(x= "Environments", y="Parameters")

stvc$CV[stvc$CV == 948.7212] <- NA
stvc$CV[stvc$CV == 130.8657] <- NA
stvc$CV[stvc$CV == 129.7948] <- NA

stvc_aft <- stvc |>  select(H2, CV, trait, env, geno_sig) |> pivot_longer(H2:CV) |>  
  ggplot(aes(x=env,y=value))+geom_col(aes(fill = trait), position = position_dodge())+ facet_grid(name~trait, scales = "free_y") + theme(legend.position = "none", axis.text.x = element_text(angle = 90)) + labs(x= "Environments", y="Parameters")

data.list = split(pheno, f = pheno$trial)

   z = data.list[["Stapleford_2018"]]
   z = data.list[["Harare_2021"]]
   table(z$`Plant lodging`)
## 
##  0  1  2  3  4  5 
## 91 38  6 22 12 11
   z = data.list[["Harare_2020"]]
   table(z$`Plant lodging`)
## 
##   0   1   2   3   4   5 
## 108  39  32  19   5  13
    z = data.list[["Harare_2022"]]
   table(z$`Plant lodging`)
## 
##  0  1  2  3  4  5 
##  9 16 10 14  6 35
   stvc |> mutate(CV = case_when(CV == 948.7212~NA))
##                                       geno        error     H2 CV         trait
## Banket_2018.Plant lodging           0.7079       1.1989 0.3713 NA Plant lodging
## Banket_2018.Plant Height          484.7230     123.1807 0.7974 NA  Plant Height
## Banket_2018.Grain Yield        165411.8384  324929.2135 0.3373 NA   Grain Yield
## Banket_2019.Plant lodging               NA           NA     NA NA Plant lodging
## Banket_2019.Plant Height                NA           NA     NA NA  Plant Height
## Banket_2019.Grain Yield        184876.1880  295453.6310 0.3849 NA   Grain Yield
## Banket_2020.Plant lodging               NA           NA     NA NA Plant lodging
## Banket_2020.Plant Height          199.5743      40.6721 0.8307 NA  Plant Height
## Banket_2020.Grain Yield        300108.5338  891371.4859 0.2519 NA   Grain Yield
## Bindura_2018.Plant lodging          1.7364       0.4628 0.7895 NA Plant lodging
## Bindura_2018.Plant Height         623.1939      88.5958 0.8755 NA  Plant Height
## Bindura_2018.Grain Yield       330922.1419  172315.2347 0.6576 NA   Grain Yield
## Bindura_2020.Plant lodging              NA           NA     NA NA Plant lodging
## Bindura_2020.Plant Height          33.9428     285.5884 0.1062 NA  Plant Height
## Bindura_2020.Grain Yield       189852.6266  697750.7882 0.2139 NA   Grain Yield
## Bindura_2021.Plant lodging              NA           NA     NA NA Plant lodging
## Bindura_2021.Plant Height               NA           NA     NA NA  Plant Height
## Bindura_2021.Grain Yield       365074.5146  250284.6771 0.5933 NA   Grain Yield
## Chisumbanje_2018.Plant lodging      1.5659       0.5192 0.7510 NA Plant lodging
## Chisumbanje_2018.Plant Height     845.3728     124.2484 0.8719 NA  Plant Height
## Chisumbanje_2018.Grain Yield   259732.3391  230272.7664 0.5301 NA   Grain Yield
## Chisumbanje_2019.Plant lodging      0.7430       0.9201 0.4468 NA Plant lodging
## Chisumbanje_2019.Plant Height     158.9274     118.1124 0.5737 NA  Plant Height
## Chisumbanje_2019.Grain Yield     9787.6001  272120.6313 0.0347 NA   Grain Yield
## Goromonzi_2018.Plant lodging            NA           NA     NA NA Plant lodging
## Goromonzi_2018.Plant Height       539.0736      67.5482 0.8886 NA  Plant Height
## Goromonzi_2018.Grain Yield     184716.1170  168639.1436 0.5227 NA   Grain Yield
## Goromonzi_2019.Plant lodging        0.2935       0.2571 0.5331 NA Plant lodging
## Goromonzi_2019.Plant Height       283.4581      90.1553 0.7587 NA  Plant Height
## Goromonzi_2019.Grain Yield     301347.8832  342784.4546 0.4678 NA   Grain Yield
## Goromonzi_2020.Plant lodging            NA           NA     NA NA Plant lodging
## Goromonzi_2020.Plant Height       140.9136      50.0625 0.7379 NA  Plant Height
## Goromonzi_2020.Grain Yield     281956.8209  414415.4873 0.4049 NA   Grain Yield
## Goromonzi_2021.Plant lodging            NA           NA     NA NA Plant lodging
## Goromonzi_2021.Plant Height       193.6497      53.6208 0.7831 NA  Plant Height
## Goromonzi_2021.Grain Yield     238728.0291   87261.9853 0.7323 NA   Grain Yield
## Goromonzi_2022.Plant lodging        0.4584       0.8481 0.3508 NA Plant lodging
## Goromonzi_2022.Plant Height       338.2353     207.7936 0.6194 NA  Plant Height
## Goromonzi_2022.Grain Yield     100029.8381  319937.0838 0.2382 NA   Grain Yield
## Goromonzi_2023.Plant lodging            NA           NA     NA NA Plant lodging
## Goromonzi_2023.Plant Height        22.6735     166.5719 0.1198 NA  Plant Height
## Goromonzi_2023.Grain Yield     538135.4862  256713.1407 0.6770 NA   Grain Yield
## Gwebi_2020.Plant lodging            1.1696       0.7089 0.6226 NA Plant lodging
## Gwebi_2020.Plant Height            97.8728     110.9833 0.4686 NA  Plant Height
## Gwebi_2020.Grain Yield          11419.4391  471465.2680 0.0236 NA   Grain Yield
## Gwebi_2021.Plant lodging            2.1307       0.2502 0.8949 NA Plant lodging
## Gwebi_2021.Plant Height           218.4241      37.0886 0.8548 NA  Plant Height
## Gwebi_2021.Grain Yield         385098.9591  468916.0707 0.4509 NA   Grain Yield
## Harare_2018.Plant lodging           0.7686       0.5908 0.5654 NA Plant lodging
## Harare_2018.Plant Height          527.0460      98.0460 0.8431 NA  Plant Height
## Harare_2018.Grain Yield        106664.3020  326694.1673 0.2461 NA   Grain Yield
## Harare_2019.Plant lodging           0.3409       1.4088 0.1948 NA Plant lodging
## Harare_2019.Plant Height           80.7317     228.7499 0.2609 NA  Plant Height
## Harare_2019.Grain Yield             0.1686 1666412.8585 0.0000 NA   Grain Yield
## Harare_2020.Plant lodging           0.0000       2.1674 0.0000 NA Plant lodging
## Harare_2020.Plant Height           99.2981     148.1606 0.4013 NA  Plant Height
## Harare_2020.Grain Yield        754984.8867  810821.7839 0.4822 NA   Grain Yield
## Harare_2021.Plant lodging           0.0586       2.5351 0.0226 NA Plant lodging
## Harare_2021.Plant Height          106.7914      30.6941 0.7767 NA  Plant Height
## Harare_2021.Grain Yield        322033.7747  460802.9537 0.4114 NA   Grain Yield
## Harare_2022.Plant lodging           1.3000       1.9318 0.4023 NA Plant lodging
## Harare_2022.Plant Height           16.6199     163.4875 0.0923 NA  Plant Height
## Harare_2022.Grain Yield         42130.5542  471560.2576 0.0820 NA   Grain Yield
## Kadoma_2018.Plant lodging           0.8203       1.1468 0.4170 NA Plant lodging
## Kadoma_2018.Plant Height           57.1665     149.5267 0.2766 NA  Plant Height
## Kadoma_2018.Grain Yield        447777.9957 1003626.5721 0.3085 NA   Grain Yield
## Kadoma_2021.Plant lodging           0.1188       1.9046 0.0587 NA Plant lodging
## Kadoma_2021.Plant Height           35.3694     311.7297 0.1019 NA  Plant Height
## Kadoma_2021.Grain Yield         33215.5894 2089801.1610 0.0156 NA   Grain Yield
## Mazowe_2022.Plant lodging           0.0152       0.4066 0.0359 NA Plant lodging
## Mazowe_2022.Plant Height          219.8637      35.4920 0.8610 NA  Plant Height
## Mazowe_2022.Grain Yield         42707.5208  161563.0444 0.2091 NA   Grain Yield
## Mthampden_2023.Plant lodging        0.2155       0.3911 0.3552 NA Plant lodging
## Mthampden_2023.Plant Height        94.3963      51.7369 0.6460 NA  Plant Height
## Mthampden_2023.Grain Yield      25746.4389  157702.2453 0.1403 NA   Grain Yield
## Mutare_2018.Plant lodging           0.4537       1.2020 0.2740 NA Plant lodging
## Mutare_2018.Plant Height          587.3138     245.3152 0.7054 NA  Plant Height
## Mutare_2018.Grain Yield        619789.7743  361342.1945 0.6317 NA   Grain Yield
## Mutare_2019.Plant lodging               NA           NA     NA NA Plant lodging
## Mutare_2019.Plant Height                NA           NA     NA NA  Plant Height
## Mutare_2019.Grain Yield        208788.9493  204126.3859 0.5056 NA   Grain Yield
## Mutare_2020.Plant lodging               NA           NA     NA NA Plant lodging
## Mutare_2020.Plant Height                NA           NA     NA NA  Plant Height
## Mutare_2020.Grain Yield        313346.1097  276250.0282 0.5315 NA   Grain Yield
## Mutare_2021.Plant lodging           0.0000       3.0276 0.0000 NA Plant lodging
## Mutare_2021.Plant Height                NA           NA     NA NA  Plant Height
## Mutare_2021.Grain Yield         99592.7483  618664.7397 0.1387 NA   Grain Yield
## Shamva_2018.Plant lodging           0.0782       0.3743 0.1727 NA Plant lodging
## Shamva_2018.Plant Height          493.6128      83.9654 0.8546 NA  Plant Height
## Shamva_2018.Grain Yield         43713.4691  150712.6572 0.2248 NA   Grain Yield
## Shamva_2019.Plant lodging           0.0954       0.4834 0.1649 NA Plant lodging
## Shamva_2019.Plant Height          174.9860    1440.9164 0.1083 NA  Plant Height
## Shamva_2019.Grain Yield         48971.9856  207604.4514 0.1909 NA   Grain Yield
## Shamva_2022.Plant lodging           0.0107       0.0218 0.3294 NA Plant lodging
## Shamva_2022.Plant Height           32.7946     114.7649 0.2222 NA  Plant Height
## Shamva_2022.Grain Yield         15230.5767  212704.1398 0.0668 NA   Grain Yield
## Stapleford_2018.Plant lodging       0.0000       0.2778 0.0000 NA Plant lodging
## Stapleford_2018.Plant Height      353.5130     232.4261 0.6033 NA  Plant Height
## Stapleford_2018.Grain Yield    649899.6569  311598.1297 0.6759 NA   Grain Yield
## Stapleford_2019.Plant lodging       0.7366       0.3767 0.6616 NA Plant lodging
## Stapleford_2019.Plant Height       81.0351      54.6182 0.5974 NA  Plant Height
## Stapleford_2019.Grain Yield    171198.8059  168382.2817 0.5041 NA   Grain Yield
## Stapleford_2020.Plant lodging           NA           NA     NA NA Plant lodging
## Stapleford_2020.Plant Height       99.5641     177.8913 0.3588 NA  Plant Height
## Stapleford_2020.Grain Yield    653844.4623  442580.0255 0.5963 NA   Grain Yield
## Stapleford_2022.Plant lodging       1.1660       0.9009 0.5641 NA Plant lodging
## Stapleford_2022.Plant Height      116.3675     113.0435 0.5072 NA  Plant Height
## Stapleford_2022.Grain Yield    311759.2963  301908.6413 0.5080 NA   Grain Yield
## Stapleford_2023.Plant lodging           NA           NA     NA NA Plant lodging
## Stapleford_2023.Plant Height      163.7645     137.8405 0.5430 NA  Plant Height
## Stapleford_2023.Grain Yield     84997.1803  319570.5776 0.2101 NA   Grain Yield
##                                             env geno_sig
## Banket_2018.Plant lodging           Banket_2018        *
## Banket_2018.Plant Height            Banket_2018        *
## Banket_2018.Grain Yield             Banket_2018        *
## Banket_2019.Plant lodging           Banket_2019     <NA>
## Banket_2019.Plant Height            Banket_2019     <NA>
## Banket_2019.Grain Yield             Banket_2019        *
## Banket_2020.Plant lodging           Banket_2020     <NA>
## Banket_2020.Plant Height            Banket_2020        *
## Banket_2020.Grain Yield             Banket_2020        *
## Bindura_2018.Plant lodging         Bindura_2018        *
## Bindura_2018.Plant Height          Bindura_2018        *
## Bindura_2018.Grain Yield           Bindura_2018        *
## Bindura_2020.Plant lodging         Bindura_2020     <NA>
## Bindura_2020.Plant Height          Bindura_2020       ns
## Bindura_2020.Grain Yield           Bindura_2020        *
## Bindura_2021.Plant lodging         Bindura_2021     <NA>
## Bindura_2021.Plant Height          Bindura_2021     <NA>
## Bindura_2021.Grain Yield           Bindura_2021        *
## Chisumbanje_2018.Plant lodging Chisumbanje_2018        *
## Chisumbanje_2018.Plant Height  Chisumbanje_2018        *
## Chisumbanje_2018.Grain Yield   Chisumbanje_2018        *
## Chisumbanje_2019.Plant lodging Chisumbanje_2019        *
## Chisumbanje_2019.Plant Height  Chisumbanje_2019        *
## Chisumbanje_2019.Grain Yield   Chisumbanje_2019       ns
## Goromonzi_2018.Plant lodging     Goromonzi_2018     <NA>
## Goromonzi_2018.Plant Height      Goromonzi_2018        *
## Goromonzi_2018.Grain Yield       Goromonzi_2018        *
## Goromonzi_2019.Plant lodging     Goromonzi_2019        *
## Goromonzi_2019.Plant Height      Goromonzi_2019        *
## Goromonzi_2019.Grain Yield       Goromonzi_2019        *
## Goromonzi_2020.Plant lodging     Goromonzi_2020     <NA>
## Goromonzi_2020.Plant Height      Goromonzi_2020        *
## Goromonzi_2020.Grain Yield       Goromonzi_2020        *
## Goromonzi_2021.Plant lodging     Goromonzi_2021     <NA>
## Goromonzi_2021.Plant Height      Goromonzi_2021        *
## Goromonzi_2021.Grain Yield       Goromonzi_2021        *
## Goromonzi_2022.Plant lodging     Goromonzi_2022        *
## Goromonzi_2022.Plant Height      Goromonzi_2022        *
## Goromonzi_2022.Grain Yield       Goromonzi_2022        *
## Goromonzi_2023.Plant lodging     Goromonzi_2023     <NA>
## Goromonzi_2023.Plant Height      Goromonzi_2023       ns
## Goromonzi_2023.Grain Yield       Goromonzi_2023        *
## Gwebi_2020.Plant lodging             Gwebi_2020        *
## Gwebi_2020.Plant Height              Gwebi_2020        *
## Gwebi_2020.Grain Yield               Gwebi_2020       ns
## Gwebi_2021.Plant lodging             Gwebi_2021        *
## Gwebi_2021.Plant Height              Gwebi_2021        *
## Gwebi_2021.Grain Yield               Gwebi_2021        *
## Harare_2018.Plant lodging           Harare_2018        *
## Harare_2018.Plant Height            Harare_2018        *
## Harare_2018.Grain Yield             Harare_2018        *
## Harare_2019.Plant lodging           Harare_2019        *
## Harare_2019.Plant Height            Harare_2019        *
## Harare_2019.Grain Yield             Harare_2019       ns
## Harare_2020.Plant lodging           Harare_2020       ns
## Harare_2020.Plant Height            Harare_2020        *
## Harare_2020.Grain Yield             Harare_2020        *
## Harare_2021.Plant lodging           Harare_2021       ns
## Harare_2021.Plant Height            Harare_2021        *
## Harare_2021.Grain Yield             Harare_2021        *
## Harare_2022.Plant lodging           Harare_2022        *
## Harare_2022.Plant Height            Harare_2022       ns
## Harare_2022.Grain Yield             Harare_2022       ns
## Kadoma_2018.Plant lodging           Kadoma_2018        *
## Kadoma_2018.Plant Height            Kadoma_2018       ns
## Kadoma_2018.Grain Yield             Kadoma_2018        *
## Kadoma_2021.Plant lodging           Kadoma_2021       ns
## Kadoma_2021.Plant Height            Kadoma_2021       ns
## Kadoma_2021.Grain Yield             Kadoma_2021       ns
## Mazowe_2022.Plant lodging           Mazowe_2022       ns
## Mazowe_2022.Plant Height            Mazowe_2022        *
## Mazowe_2022.Grain Yield             Mazowe_2022       ns
## Mthampden_2023.Plant lodging     Mthampden_2023        *
## Mthampden_2023.Plant Height      Mthampden_2023        *
## Mthampden_2023.Grain Yield       Mthampden_2023       ns
## Mutare_2018.Plant lodging           Mutare_2018        *
## Mutare_2018.Plant Height            Mutare_2018        *
## Mutare_2018.Grain Yield             Mutare_2018        *
## Mutare_2019.Plant lodging           Mutare_2019     <NA>
## Mutare_2019.Plant Height            Mutare_2019     <NA>
## Mutare_2019.Grain Yield             Mutare_2019        *
## Mutare_2020.Plant lodging           Mutare_2020     <NA>
## Mutare_2020.Plant Height            Mutare_2020     <NA>
## Mutare_2020.Grain Yield             Mutare_2020        *
## Mutare_2021.Plant lodging           Mutare_2021       ns
## Mutare_2021.Plant Height            Mutare_2021     <NA>
## Mutare_2021.Grain Yield             Mutare_2021       ns
## Shamva_2018.Plant lodging           Shamva_2018       ns
## Shamva_2018.Plant Height            Shamva_2018        *
## Shamva_2018.Grain Yield             Shamva_2018        *
## Shamva_2019.Plant lodging           Shamva_2019        *
## Shamva_2019.Plant Height            Shamva_2019       ns
## Shamva_2019.Grain Yield             Shamva_2019        *
## Shamva_2022.Plant lodging           Shamva_2022        *
## Shamva_2022.Plant Height            Shamva_2022       ns
## Shamva_2022.Grain Yield             Shamva_2022       ns
## Stapleford_2018.Plant lodging   Stapleford_2018       ns
## Stapleford_2018.Plant Height    Stapleford_2018        *
## Stapleford_2018.Grain Yield     Stapleford_2018        *
## Stapleford_2019.Plant lodging   Stapleford_2019        *
## Stapleford_2019.Plant Height    Stapleford_2019        *
## Stapleford_2019.Grain Yield     Stapleford_2019        *
## Stapleford_2020.Plant lodging   Stapleford_2020     <NA>
## Stapleford_2020.Plant Height    Stapleford_2020        *
## Stapleford_2020.Grain Yield     Stapleford_2020        *
## Stapleford_2022.Plant lodging   Stapleford_2022        *
## Stapleford_2022.Plant Height    Stapleford_2022        *
## Stapleford_2022.Grain Yield     Stapleford_2022        *
## Stapleford_2023.Plant lodging   Stapleford_2023     <NA>
## Stapleford_2023.Plant Height    Stapleford_2023        *
## Stapleford_2023.Grain Yield     Stapleford_2023        *

Bayesian model

# Filter env which has significance to genotype effect
stmt = readRDS(file = "Saves/stmt.rds")

diagnostic = do.call(rbind, lapply(stmt, function(x){
  aa = as.data.frame(x$summary) |> rownames_to_column("trait")
  aa
})) |> rownames_to_column("env") |> 
  mutate_at("env", str_replace, "\\..*", "") |> 
  dplyr::select(env, trait, geno) |> 
  separate(geno, into = c("genovar", "sig"), sep = " ")

# diagnostic
# dim(pheno)
# dim(pheno_filt)
# unique(pheno_filt$trial)
# levels(pheno_filt$trial)


# inpheno# install.packages("ProbBreed")
library(ProbBreed)
# head(pheno_filt)

i=traits[3]

# Getting the path of your current open file
# current_path = rstudioapi::getActiveDocumentContext()$path 
# setwd(dirname(current_path ))
# print( getwd() )
# head(pheno)
for (i in traits) {
  message(paste("Analysis to", i))
  # print(paste("Analysis to", i)) 
  # cat("====> Trait:", i, fill = TRUE)
  pheno_filt <- droplevels(pheno[
    which(pheno$trial %in% diagnostic[which(diagnostic$trait == i &
                                               diagnostic$sig == "*"), "env"]), ])
  

   
  dim(pheno_filt)
  mod = bayes_met(data = pheno_filt,
                gen = "geno",
                loc = "LOCATION_NAME",
                repl = c("rep"),
                trait = paste(i),
                reg = NULL,
                year = "CROP_SEASON",
                res.het = TRUE,
                iter = 4000, cores = 4, chain = 4)
  
#save(mod,file = paste0("Saves/",i,".rda"))  
saveRDS(mod,file = paste0("Saves/",i,"_year.rds"))
rm(mod)
}

rm(list=ls())

mod_PH=readRDS("Saves/Plant Height.rds")
mod_PL=readRDS("Saves/Plant lodging.rds")
mod_GY=readRDS("Saves/Grain Yield.rds")

 

# setwd("./Saves/PH")
res_PH = prob_sup(extr = outs_PH,
                   int = .2,
                   increase = FALSE,
                   save.df = FALSE,
                   verbose = FALSE)
# setwd("../PL")

res_PL = prob_sup(extr = outs_PL,
                   int = .2,
                   increase = FALSE,
                   save.df = FALSE,
                   verbose = FALSE)
# setwd("../GY")

res_GY = prob_sup(extr = outs_GY,
                   int = .2,
                   increase = TRUE,
                   save.df = FALSE,
                   verbose = FALSE)

# setwd("../")

save(res_PL,file = "Saves/res_PL_year.rda")
save(res_PH,file = "Saves/res_PH_year.rda")
save(res_GY,file = "Saves/res_GY_year.rda")

Variances

outs_PH=readRDS("Saves/outs_PH.rds")
outs_PL=readRDS("Saves/outs_PL.rds")
outs_GY=readRDS("Saves/outs_GY.rds")


outs_PH$variances$trait="PH"
outs_PL$variances$trait="PL"
outs_GY$variances$trait="GY"

varcomp=rbind(outs_GY$variances,outs_PH$variances,outs_PL$variances) 

library(tidyverse)
head(varcomp)
varcomp |> ggplot(aes(x=effect,y=var, fill=trait)) + geom_col(position = "dodge") + theme(axis.text.x = element_text(angle=90), legend.position = "none") + facet_wrap(~trait,scales="free")
outs_PH=readRDS("Saves/outs_PH_year.rds")
outs_PL=readRDS("Saves/outs_PL_year.rds")
outs_GY=readRDS("Saves/outs_GY_year.rds")


outs_PH$variances$trait="Plant Height"
outs_PL$variances$trait="Plant Lodging"
outs_GY$variances$trait="Grain Yield"

varcomp=rbind(outs_GY$variances,outs_PH$variances,outs_PL$variances) 

library(tidyverse)
head(varcomp)
vc=varcomp |> ggplot(aes(x=effect,y=var, fill=trait)) + geom_col(position = "dodge") + theme(axis.text.x = element_text(angle=90), legend.position = "none") +
 geom_errorbar(aes(ymin = var-sd,ymax=var+sd),width=0.1)+ facet_wrap(~trait,scales="free")

ggsave(plot=vc,device = "pdf",filename ="Plots/varcomp1.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=vc,device = "png",filename ="Plots/varcomp1.png",dpi = "retina",height = 6,width = 8,bg = "white")

Density

library(ProbBreed)
library(tidyverse)
outs_PH=readRDS("Saves/outs_PH_year.rds")
outs_PL=readRDS("Saves/outs_PL_year.rds")
outs_GY=readRDS("Saves/outs_GY_year.rds")


outs_PH$ppcheck
outs_PL$ppcheck
outs_GY$ppcheck

# plot(outs_PH, category = "density")
# plot(outs_PH, category = "histogram")

# plot(outs_PL, category = "density")
# plot(outs_PL, category = "histogram")

# plot(outs_GY, category = "density")
# plot(outs_GY, category = "histogram")
# Density: Empirical vs Sampled phenotypic values
PL=plot(outs_PL)+ggtitle("")+ labs(x="Plant Lodging")
PH=plot(outs_PH)+ggtitle("")
GY=plot(outs_GY)+ggtitle("")
library(ggpubr)
plot.density = ggarrange(GY,PH,PL,labels = c("A","B","C"),nrow = 3, common.legend = T)


ggsave(plot=plot.density,device = "pdf",filename ="Plots/plots_density.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=plot.density,device = "png",filename ="Plots/plots_density.png",dpi = "retina",height = 6,width = 8,bg = "white")

Probsup within locations

rm(list=ls())


library(ProbBreed)
library(tidyverse)
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")



df=read.csv("BPSI/bpsi.csv")

gen.sel = df[which(df$sel=="selected"),"gen"]

GY=plot(res_GY, category = "perfo", level = "within")
PH=plot(res_PH, category = "perfo", level = "within")
PL=plot(res_PL, category = "perfo", level = "within")

unique(GY$data$fac)
cols.label=c("Location", "Year")
names(cols.label)=c(unique(GY$data$fac)) 




labels <- function(x) {
  gsub("G_", "V0", x) 
}



GY=GY+ theme(axis.text.y = element_text(size=4,color = ifelse(GY$data$gen %in% gen.sel,  "#E41A1C", "black"))) + labs( title = "Grain Yield")+facet_wrap(~fac,scales="free", labeller = labeller(.cols=cols.label))+ scale_y_discrete(labels = labels)

PH=PH+ theme(axis.text.y = element_text(size=4, color = ifelse(PH$data$gen %in% gen.sel,  "#E41A1C", "black"))) + labs( title = "Plant Height")+facet_wrap(~fac,scales="free", labeller = labeller(.cols=cols.label))+ scale_y_discrete(labels = labels)
PL=PL+ theme(axis.text.y = element_text(size=4, color = ifelse(PL$data$gen %in% gen.sel,  "#E41A1C", "black"))) + labs( title = "Plant Lodging")+facet_wrap(~fac,scales="free", labeller = labeller(.cols=cols.label))+ scale_y_discrete(labels = labels)

library(ggpubr)
plot.within = ggarrange(GY,PH,PL,labels = c("A","B","C"),ncol=3, common.legend = T)


ggsave(plot=plot.within,device = "pdf",filename ="Plots/plots_perfo_within.pdf",dpi = "retina",height = 8,width = 10)
ggsave(plot=plot.within,device = "png",filename ="Plots/plots_perfo_within.png",dpi = "retina",height = 8,width = 10,bg = "white")

Probsup across locations

library(ProbBreed)
library(tidyverse)
load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")



df=read.csv("BPSI/bpsi.csv")

gen.sel = df[which(df$sel=="selected"),"gen"]


PH=plot(res_PH)
PL=plot(res_PL)
GY=plot(res_GY)

labels <- function(x) {
  gsub("G_", "V0", x) 
}


library(RColorBrewer)
RColorBrewer::brewer.pal(n=8,name = "Spectral")
display.brewer.all(n=8)

PH=PH +   geom_point(aes(  x = PH$data$ID,y = PH$data$prob,
fill = ifelse(PH$data$ID %in% gen.sel,  "#377EB8","#E41A1C" )) 
,size = 2,shape = 21,stroke = 0.4) +guides(fill = "none") +  
  labs(y = "", x = "", title = "Plant Height")+theme(axis.text.x = element_text(size = 5, angle = 90, hjust = 1, vjust = 0.5))+ scale_x_discrete(labels = labels)

PL=PL +   geom_point(aes(  x = PL$data$ID,y = PL$data$prob,
fill = ifelse(PL$data$ID %in% gen.sel,  "#377EB8","#E41A1C" )) 
,size = 2,shape = 21,stroke = 0.4) +guides(fill = "none") +  
  labs(y = "", x = "", title = "Plant Lodging")+theme(axis.text.x = element_text(size = 5, angle = 90, hjust = 1, vjust = 0.5))+ scale_x_discrete(labels = labels)

GY=GY +   geom_point(aes(  x = GY$data$ID,y = GY$data$prob,
fill = ifelse(GY$data$ID %in% gen.sel,  "#377EB8","#E41A1C" )) 
,size = 2,shape = 21,stroke = 0.4) +guides(fill = "none") +  
  labs(y = "", x = "", title = "Grain Yield")+theme(axis.text.x = element_text(size = 5, angle = 90, hjust = 1, vjust = 0.5))+ scale_x_discrete(labels = labels)


library(ggpubr)
plot.perfo = ggarrange(GY,PH,PL,labels = c("A","B","C"),nrow = 3, common.legend = T)


ggsave(plot=plot.perfo,device = "pdf",filename ="Plots/plots_perfo.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=plot.perfo,device = "png",filename ="Plots/plots_perfo.png",dpi = "retina",height = 6,width = 8)

## Probsup across locations pairwise

load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")

df=read.csv("BPSI/bpsi.csv")

gen.sel = df[which(df$sel=="selected"),"gen"]


GY=plot(res_GY, category = "pair_perfo")
PH=plot(res_PH, category = "pair_perfo")
PL=plot(res_PL, category = "pair_perfo")

library(tidyverse)



target2 = gen.sel
# load("margsarq.rda")
# load("b.rda")
# load("a.rda")
a=GY
b=PH
c=PL


a$data=a$data[which(a$data$x  %in% target2   ),]
b$data=b$data[which(b$data$x  %in% target2   ),]
c$data=c$data[which(c$data$x  %in% target2  ),]
a$data$model ="GY"
b$data$model ="PH"
c$data$model = "PL"
dim(a$data)

highlighted_obs <- c(gen.sel)



GY=a +
  theme(
          panel.background = element_blank(), legend.position = "top", legend.direction = "horizontal",
          axis.text.y = element_text(color = ifelse(  a$data$y %in% highlighted_obs , "red", "black"), size = 5), strip.text = element_text(size = 16)) +
          scale_fill_viridis_c(direction = 1, na.value = "white",
            limits = c(0, 1), option = "turbo") + guides(fill = guide_colorbar(barwidth = 7, barheight = 1.5, title.position = "top", title.hjust = 0.5)) + labs(title = "Grain Yield",y="")+ scale_x_discrete(labels = labels)+ scale_y_discrete(labels = labels)

# b$data=b$data[which(b$data$x  %in% target2 &  b$data$y  %in% target2),]

PH=b +
  theme(
          panel.background = element_blank(), legend.position = "top", legend.direction = "horizontal",
          axis.text.y = element_text(color = ifelse(  b$data$y %in% highlighted_obs , "red", "black"), size = 5), strip.text = element_text(size = 16)) +
          scale_fill_viridis_c(direction = 1, na.value = "white",
            limits = c(0, 1), option = "turbo") + guides(fill = guide_colorbar(barwidth = 7, barheight = 1.5, title.position = "top", title.hjust = 0.5)) + labs(title = "Plant Heigth", x="")+ scale_x_discrete(labels = labels)+ scale_y_discrete(labels = labels)

PL=c +
  theme(    panel.background = element_blank(), legend.position = "top", legend.direction = "horizontal",
          axis.text.y = element_text(color = ifelse(  c$data$y %in% highlighted_obs , "red", "black"), size = 5), strip.text = element_text(size = 16)) +
          scale_fill_viridis_c(direction = 1, na.value = "white",
            limits = c(0, 1), option = "turbo") + guides(fill = guide_colorbar(barwidth = 7, barheight = 1.5, title.position = "top", title.hjust = 0.5)) + labs(title = "Plant Lodging",x="",y="")+ scale_x_discrete(labels = labels)+ scale_y_discrete(labels = labels)




plot.pair = ggarrange(GY,PH,PL,labels = c("A","B","C"),ncol = 3, common.legend = T)






plot.pair = ggarrange(GY,PH,PL,labels = c("A","B","C"),ncol = 3, common.legend = T)


ggsave(plot=plot.pair,device = "pdf",filename ="Plots/plots_pair.pdf",dpi = "retina",height = 6,width = 8)
ggsave(plot=plot.pair,device = "png",filename ="Plots/plots_pair.png",dpi = "retina",height = 6,width = 8, bg = "white")

# BPSI - w 2,1,1

load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")

source("data/bpsi_fun.R")

models= vector("list",length(3))


models[[1]] = res_GY;
models[[2]] = res_PH
models[[3]] = res_PL
models[[3]]$across$perfo <- models[[3]]$across$perfo[-which(models[[3]]$across$perfo$ID=="G_55"),]
names(models) <- c("GY","PH","PL")
BPSI_soy=BPSI(modlist = models,increase =c(TRUE,FALSE,FALSE),omega = c(2,1,1),int = 0.1,save.df = T,verbose = T )
## -> Genotypes selected using BPSI
## Warning in dir.create(path = paste0(getwd(), "/BPSI")):
## '/Users/tiagobchagas/Documents/UFV/Zimbabwe_trials/BPSI' already exists
## Process completed!
df=print(BPSI_soy)
## ==> Considering an intensity of 10%, here are the selected candidates:
##      GY   PH   PL  bpsi      sel  gen
## 1   0.5  4.0  4.0   8.5 selected G_20
## 2   8.0  6.0  4.0  18.0 selected  G_3
## 3   3.0  2.0 17.5  22.5 selected G_38
## 4   7.0  5.0 18.0  30.0 selected G_45
## 5  13.0  3.0 15.0  31.0 selected  G_4
## 6   7.0  8.0 24.5  39.5 selected G_13
## 7  29.0 12.0  5.0  46.0 selected G_16
## 8  19.0 18.0 14.0  51.0 selected  G_1
## 9   2.0 49.0  1.5  52.5 selected G_32
## 10 19.0 21.5 12.0  52.5 selected G_33
## 11 27.0  7.0 19.0  53.0  not_sel G_14
## 12 15.0 38.0  4.5  57.5  not_sel G_21
## 13 10.0 17.0 33.0  60.0  not_sel G_31
## 14 57.0  1.0  2.5  60.5  not_sel G_54
## 15  2.0 52.0  7.0  61.0  not_sel G_50
## 16 17.0 17.5 28.0  62.5  not_sel G_28
## 17 48.0  3.5 11.0  62.5  not_sel G_64
## 18  5.0 47.0 13.0  65.0  not_sel G_93
## 19 13.0 42.0 13.5  68.5  not_sel G_40
## 20 34.0 24.0 11.5  69.5  not_sel G_49
## 21 23.0 25.0 23.0  71.0  not_sel G_60
## 22 37.0 18.0 17.0  72.0  not_sel G_39
## 23 20.5 39.0 13.0  72.5  not_sel G_68
## 24 36.0 29.0  8.0  73.0  not_sel  G_8
## 25 24.0 34.0 18.5  76.5  not_sel G_77
## 26  3.0 61.0 20.0  84.0  not_sel G_53
## 27 18.0 24.0 47.0  89.0  not_sel G_36
## 28 21.0 38.0 30.0  89.0  not_sel G_47
## 29 33.0 20.5 38.0  91.5  not_sel G_52
## 30 23.0 30.0 39.0  92.0  not_sel G_29
## 31  9.0 28.0 55.0  92.0  not_sel G_44
## 32 14.0 76.0  2.0  92.0  not_sel  G_9
## 33 32.0  9.5 51.0  92.5  not_sel G_11
## 34 43.0 28.0 21.5  92.5  not_sel G_74
## 35 11.0 58.0 25.0  94.0  not_sel G_42
## 36 44.0 13.0 43.5 100.5  not_sel G_58
## 37 12.0 36.0 53.0 101.0  not_sel G_25
## 38 73.0 10.0 24.0 107.0  not_sel G_67
## 39 11.0 76.0 20.5 107.5  not_sel G_46
## 40 68.0 23.0 18.0 109.0  not_sel G_10
## 41 55.0  4.5 50.0 109.5  not_sel  G_2
## 42 49.0 16.0 45.0 110.0  not_sel G_30
## 43 77.0 13.0 21.0 111.0  not_sel G_75
## 44 15.0 53.0 44.0 112.0  not_sel G_92
## 45 54.0  5.0 54.0 113.0  not_sel G_17
## 46 10.0 76.0 29.5 115.5  not_sel G_43
## 47 81.0  5.5 29.0 115.5  not_sel  G_5
## 48 31.0 25.0 60.0 116.0  not_sel G_89
## 49 33.5 22.0 63.0 118.5  not_sel  G_6
## 50 65.0 15.0 40.0 120.0  not_sel G_88
## 51  8.0 76.0 40.0 124.0  not_sel G_34
## 52 30.5 72.0 22.0 124.5  not_sel G_70
## 53 35.0 76.0 16.0 127.0  not_sel G_63
## 54 51.0 76.0  0.5 127.5  not_sel G_27
## 55 47.0 15.5 66.0 128.5  not_sel G_94
## 56 45.0 60.0 24.0 129.0  not_sel G_35
## 57 50.0 76.0  3.0 129.0  not_sel G_51
## 58 20.0 76.0 34.0 130.0  not_sel G_26
## 59 59.0 22.5 52.0 133.5  not_sel G_61
## 60 42.0 38.0 57.0 137.0  not_sel G_91
## 61 39.0 38.0 64.0 141.0  not_sel G_41
## 62 86.0 16.0 41.0 143.0  not_sel G_19
## 63 31.0 40.0 76.0 147.0  not_sel G_84
## 64 80.0 27.0 42.0 149.0  not_sel G_97
## 65 32.0 27.0 90.0 149.0  not_sel G_98
## 66 63.0 23.0 65.0 151.0  not_sel G_72
## 67 38.5 21.0 93.0 152.5  not_sel G_65
## 68 33.0 59.0 62.0 154.0  not_sel G_81
## 69 56.0 64.0 35.0 155.0  not_sel G_85
## 70 12.5 76.0 67.0 155.5  not_sel G_73
## 71 58.0 64.0 39.5 161.5  not_sel G_82
## 72 60.0 64.0 41.0 165.0  not_sel G_90
## 73 35.5 72.0 58.0 165.5  not_sel G_48
## 74 40.5 33.0 93.0 166.5  not_sel G_57
## 75 43.0 69.0 56.0 168.0  not_sel G_87
## 76 26.0 55.0 88.0 169.0  not_sel G_76
## 77 53.0 28.5 89.0 170.5  not_sel G_80
## 78 42.0 64.0 68.0 174.0  not_sel G_12
## 79 34.5 62.0 84.0 180.5  not_sel G_95
## 80 90.0 76.0 15.5 181.5  not_sel G_24
## 81 37.5 70.0 75.0 182.5  not_sel G_18
## 82 84.0 38.0 61.0 183.0  not_sel G_22
## 83 86.0 63.0 35.5 184.5  not_sel G_71
## 84 86.0 25.5 73.0 184.5  not_sel G_83
## 85 70.0 38.0 77.0 185.0  not_sel G_78
## 86 79.0 18.5 90.0 187.5  not_sel G_56
## 87 74.0 70.0 45.0 189.0  not_sel G_96
## 88 81.0 76.0 36.0 193.0  not_sel G_66
## 89 90.0 22.0 81.0 193.0  not_sel  G_7
## 90 72.0 36.0 86.0 194.0  not_sel G_86
## 91 45.0 76.0 74.0 195.0  not_sel G_23
## 92 90.0 38.0 69.0 197.0  not_sel G_59
## 93 45.0 76.0 78.0 199.0  not_sel G_37
## 94 38.0 76.0 85.0 199.0  not_sel G_79
## 95 45.0 64.0 93.0 202.0  not_sel G_15
## 96 90.0 76.0 46.5 212.5  not_sel G_69
## 97 45.0 76.0 93.0 214.0  not_sel G_62
gen.sel = df[which(df$sel=="selected"),"gen"]

real.names=distinct(pheno[,c(5,15)])

df <- left_join(df, real.names, by = c("gen" = "geno"))

BPSI plot - w 2,1,1

Highlighting the ranking of probability of superior performance per trait of the selected families.

library(knitr)

plot(BPSI_soy,category = "Ranks")
Ranking of probability of superior performance per trait - w (2,1,1)

Ranking of probability of superior performance per trait - w (2,1,1)

plot(BPSI_soy,category = "BPSI")
## Warning: Could not calculate the predicate for layer 3; ignored
Selected families based on BPSI rank - w (2,1,1)

Selected families based on BPSI rank - w (2,1,1)

library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
BPSI_soy$BPSI |> kable(format = "html", booktabs = TRUE, longtable=TRUE)
GY PH PL bpsi sel gen
G_20 0.5 4.0 4.0 8.5 selected G_20
G_3 8.0 6.0 4.0 18.0 selected G_3
G_38 3.0 2.0 17.5 22.5 selected G_38
G_45 7.0 5.0 18.0 30.0 selected G_45
G_4 13.0 3.0 15.0 31.0 selected G_4
G_13 7.0 8.0 24.5 39.5 selected G_13
G_16 29.0 12.0 5.0 46.0 selected G_16
G_1 19.0 18.0 14.0 51.0 selected G_1
G_32 2.0 49.0 1.5 52.5 selected G_32
G_33 19.0 21.5 12.0 52.5 selected G_33
G_14 27.0 7.0 19.0 53.0 not_sel G_14
G_21 15.0 38.0 4.5 57.5 not_sel G_21
G_31 10.0 17.0 33.0 60.0 not_sel G_31
G_54 57.0 1.0 2.5 60.5 not_sel G_54
G_50 2.0 52.0 7.0 61.0 not_sel G_50
G_28 17.0 17.5 28.0 62.5 not_sel G_28
G_64 48.0 3.5 11.0 62.5 not_sel G_64
G_93 5.0 47.0 13.0 65.0 not_sel G_93
G_40 13.0 42.0 13.5 68.5 not_sel G_40
G_49 34.0 24.0 11.5 69.5 not_sel G_49
G_60 23.0 25.0 23.0 71.0 not_sel G_60
G_39 37.0 18.0 17.0 72.0 not_sel G_39
G_68 20.5 39.0 13.0 72.5 not_sel G_68
G_8 36.0 29.0 8.0 73.0 not_sel G_8
G_77 24.0 34.0 18.5 76.5 not_sel G_77
G_53 3.0 61.0 20.0 84.0 not_sel G_53
G_36 18.0 24.0 47.0 89.0 not_sel G_36
G_47 21.0 38.0 30.0 89.0 not_sel G_47
G_52 33.0 20.5 38.0 91.5 not_sel G_52
G_29 23.0 30.0 39.0 92.0 not_sel G_29
G_44 9.0 28.0 55.0 92.0 not_sel G_44
G_9 14.0 76.0 2.0 92.0 not_sel G_9
G_11 32.0 9.5 51.0 92.5 not_sel G_11
G_74 43.0 28.0 21.5 92.5 not_sel G_74
G_42 11.0 58.0 25.0 94.0 not_sel G_42
G_58 44.0 13.0 43.5 100.5 not_sel G_58
G_25 12.0 36.0 53.0 101.0 not_sel G_25
G_67 73.0 10.0 24.0 107.0 not_sel G_67
G_46 11.0 76.0 20.5 107.5 not_sel G_46
G_10 68.0 23.0 18.0 109.0 not_sel G_10
G_2 55.0 4.5 50.0 109.5 not_sel G_2
G_30 49.0 16.0 45.0 110.0 not_sel G_30
G_75 77.0 13.0 21.0 111.0 not_sel G_75
G_92 15.0 53.0 44.0 112.0 not_sel G_92
G_17 54.0 5.0 54.0 113.0 not_sel G_17
G_43 10.0 76.0 29.5 115.5 not_sel G_43
G_5 81.0 5.5 29.0 115.5 not_sel G_5
G_89 31.0 25.0 60.0 116.0 not_sel G_89
G_6 33.5 22.0 63.0 118.5 not_sel G_6
G_88 65.0 15.0 40.0 120.0 not_sel G_88
G_34 8.0 76.0 40.0 124.0 not_sel G_34
G_70 30.5 72.0 22.0 124.5 not_sel G_70
G_63 35.0 76.0 16.0 127.0 not_sel G_63
G_27 51.0 76.0 0.5 127.5 not_sel G_27
G_94 47.0 15.5 66.0 128.5 not_sel G_94
G_35 45.0 60.0 24.0 129.0 not_sel G_35
G_51 50.0 76.0 3.0 129.0 not_sel G_51
G_26 20.0 76.0 34.0 130.0 not_sel G_26
G_61 59.0 22.5 52.0 133.5 not_sel G_61
G_91 42.0 38.0 57.0 137.0 not_sel G_91
G_41 39.0 38.0 64.0 141.0 not_sel G_41
G_19 86.0 16.0 41.0 143.0 not_sel G_19
G_84 31.0 40.0 76.0 147.0 not_sel G_84
G_97 80.0 27.0 42.0 149.0 not_sel G_97
G_98 32.0 27.0 90.0 149.0 not_sel G_98
G_72 63.0 23.0 65.0 151.0 not_sel G_72
G_65 38.5 21.0 93.0 152.5 not_sel G_65
G_81 33.0 59.0 62.0 154.0 not_sel G_81
G_85 56.0 64.0 35.0 155.0 not_sel G_85
G_73 12.5 76.0 67.0 155.5 not_sel G_73
G_82 58.0 64.0 39.5 161.5 not_sel G_82
G_90 60.0 64.0 41.0 165.0 not_sel G_90
G_48 35.5 72.0 58.0 165.5 not_sel G_48
G_57 40.5 33.0 93.0 166.5 not_sel G_57
G_87 43.0 69.0 56.0 168.0 not_sel G_87
G_76 26.0 55.0 88.0 169.0 not_sel G_76
G_80 53.0 28.5 89.0 170.5 not_sel G_80
G_12 42.0 64.0 68.0 174.0 not_sel G_12
G_95 34.5 62.0 84.0 180.5 not_sel G_95
G_24 90.0 76.0 15.5 181.5 not_sel G_24
G_18 37.5 70.0 75.0 182.5 not_sel G_18
G_22 84.0 38.0 61.0 183.0 not_sel G_22
G_71 86.0 63.0 35.5 184.5 not_sel G_71
G_83 86.0 25.5 73.0 184.5 not_sel G_83
G_78 70.0 38.0 77.0 185.0 not_sel G_78
G_56 79.0 18.5 90.0 187.5 not_sel G_56
G_96 74.0 70.0 45.0 189.0 not_sel G_96
G_66 81.0 76.0 36.0 193.0 not_sel G_66
G_7 90.0 22.0 81.0 193.0 not_sel G_7
G_86 72.0 36.0 86.0 194.0 not_sel G_86
G_23 45.0 76.0 74.0 195.0 not_sel G_23
G_59 90.0 38.0 69.0 197.0 not_sel G_59
G_37 45.0 76.0 78.0 199.0 not_sel G_37
G_79 38.0 76.0 85.0 199.0 not_sel G_79
G_15 45.0 64.0 93.0 202.0 not_sel G_15
G_69 90.0 76.0 46.5 212.5 not_sel G_69
G_62 45.0 76.0 93.0 214.0 not_sel G_62

BPSI - w 3,1,1

load("Saves/res_PL_year.rda")
load("Saves/res_PH_year.rda")
load("Saves/res_GY_year.rda")

source("data/bpsi_fun.R")

models= vector("list",length(3))


models[[1]] = res_GY;
models[[2]] = res_PH
models[[3]] = res_PL
models[[3]]$across$perfo <- models[[3]]$across$perfo[-which(models[[3]]$across$perfo$ID=="G_55"),]
names(models) <- c("GY","PH","PL")
BPSI_soy=BPSI(modlist = models,increase =c(TRUE,FALSE,FALSE),omega = c(3,1,1),int = 0.1,save.df = T,verbose = T )
## -> Genotypes selected using BPSI
## Warning in dir.create(path = paste0(getwd(), "/BPSI")):
## '/Users/tiagobchagas/Documents/UFV/Zimbabwe_trials/BPSI' already exists
## Process completed!
df=print(BPSI_soy)
## ==> Considering an intensity of 10%, here are the selected candidates:
##            GY        PH         PL       bpsi      sel  gen
## 1   0.3333333  4.000000  4.0000000   8.333333 selected G_20
## 2   3.0000000  2.000000 11.6666667  16.666667 selected G_38
## 3   8.0000000  6.000000  2.6666667  16.666667 selected  G_3
## 4   8.6666667  3.000000 15.0000000  26.666667 selected  G_4
## 5   4.6666667  5.000000 18.0000000  27.666667 selected G_45
## 6   7.0000000  8.000000 16.3333333  31.333333 selected G_13
## 7  29.0000000 12.000000  3.3333333  44.333333 selected G_16
## 8  12.6666667 18.000000 14.0000000  44.666667 selected  G_1
## 9  19.0000000 14.333333 12.0000000  45.333333 selected G_33
## 10 27.0000000  4.666667 19.0000000  50.666667 selected G_14
## 11  2.0000000 49.000000  1.0000000  52.000000  not_sel G_32
## 12 15.0000000 38.000000  3.0000000  56.000000  not_sel G_21
## 13 17.0000000 11.666667 28.0000000  56.666667  not_sel G_28
## 14  6.6666667 17.000000 33.0000000  56.666667  not_sel G_31
## 15 57.0000000  1.000000  1.6666667  59.666667  not_sel G_54
## 16  1.3333333 52.000000  7.0000000  60.333333  not_sel G_50
## 17  5.0000000 47.000000  8.6666667  60.666667  not_sel G_93
## 18 48.0000000  2.333333 11.0000000  61.333333  not_sel G_64
## 19 23.0000000 25.000000 15.3333333  63.333333  not_sel G_60
## 20 13.0000000 42.000000  9.0000000  64.000000  not_sel G_40
## 21 13.6666667 39.000000 13.0000000  65.666667  not_sel G_68
## 22 34.0000000 24.000000  7.6666667  65.666667  not_sel G_49
## 23 37.0000000 12.000000 17.0000000  66.000000  not_sel G_39
## 24 24.0000000 34.000000 12.3333333  70.333333  not_sel G_77
## 25 36.0000000 29.000000  5.3333333  70.333333  not_sel  G_8
## 26 21.0000000 25.333333 30.0000000  76.333333  not_sel G_47
## 27 18.0000000 16.000000 47.0000000  81.000000  not_sel G_36
## 28  9.0000000 18.666667 55.0000000  82.666667  not_sel G_44
## 29  2.0000000 61.000000 20.0000000  83.000000  not_sel G_53
## 30 15.3333333 30.000000 39.0000000  84.333333  not_sel G_29
## 31 33.0000000 13.666667 38.0000000  84.666667  not_sel G_52
## 32 43.0000000 28.000000 14.3333333  85.333333  not_sel G_74
## 33 44.0000000 13.000000 29.0000000  86.000000  not_sel G_58
## 34  9.3333333 76.000000  2.0000000  87.333333  not_sel  G_9
## 35 12.0000000 24.000000 53.0000000  89.000000  not_sel G_25
## 36 32.0000000  6.333333 51.0000000  89.333333  not_sel G_11
## 37  7.3333333 58.000000 25.0000000  90.333333  not_sel G_42
## 38 11.0000000 76.000000 13.6666667 100.666667  not_sel G_46
## 39 68.0000000 23.000000 12.0000000 103.000000  not_sel G_10
## 40 73.0000000  6.666667 24.0000000 103.666667  not_sel G_67
## 41 49.0000000 10.666667 45.0000000 104.666667  not_sel G_30
## 42 10.0000000 76.000000 19.6666667 105.666667  not_sel G_43
## 43 77.0000000  8.666667 21.0000000 106.666667  not_sel G_75
## 44 65.0000000 15.000000 26.6666667 106.666667  not_sel G_88
## 45 10.0000000 53.000000 44.0000000 107.000000  not_sel G_92
## 46 22.3333333 22.000000 63.0000000 107.333333  not_sel  G_6
## 47 31.0000000 16.666667 60.0000000 107.666667  not_sel G_89
## 48 55.0000000  3.000000 50.0000000 108.000000  not_sel  G_2
## 49 54.0000000  3.333333 54.0000000 111.333333  not_sel G_17
## 50 81.0000000  3.666667 29.0000000 113.666667  not_sel  G_5
## 51 20.3333333 72.000000 22.0000000 114.333333  not_sel G_70
## 52 45.0000000 60.000000 16.0000000 121.000000  not_sel G_35
## 53  5.3333333 76.000000 40.0000000 121.333333  not_sel G_34
## 54 35.0000000 76.000000 10.6666667 121.666667  not_sel G_63
## 55 13.3333333 76.000000 34.0000000 123.333333  not_sel G_26
## 56 47.0000000 10.333333 66.0000000 123.333333  not_sel G_94
## 57 42.0000000 25.333333 57.0000000 124.333333  not_sel G_91
## 58 59.0000000 15.000000 52.0000000 126.000000  not_sel G_61
## 59 51.0000000 76.000000  0.3333333 127.333333  not_sel G_27
## 60 50.0000000 76.000000  2.0000000 128.000000  not_sel G_51
## 61 39.0000000 25.333333 64.0000000 128.333333  not_sel G_41
## 62 86.0000000 16.000000 27.3333333 129.333333  not_sel G_19
## 63 20.6666667 40.000000 76.0000000 136.666667  not_sel G_84
## 64 21.3333333 27.000000 90.0000000 138.333333  not_sel G_98
## 65 25.6666667 21.000000 93.0000000 139.666667  not_sel G_65
## 66 80.0000000 18.000000 42.0000000 140.000000  not_sel G_97
## 67 22.0000000 59.000000 62.0000000 143.000000  not_sel G_81
## 68 63.0000000 15.333333 65.0000000 143.333333  not_sel G_72
## 69 56.0000000 64.000000 23.3333333 143.333333  not_sel G_85
## 70 58.0000000 64.000000 26.3333333 148.333333  not_sel G_82
## 71  8.3333333 76.000000 67.0000000 151.333333  not_sel G_73
## 72 60.0000000 64.000000 27.3333333 151.333333  not_sel G_90
## 73 27.0000000 33.000000 93.0000000 153.000000  not_sel G_57
## 74 23.6666667 72.000000 58.0000000 153.666667  not_sel G_48
## 75 28.6666667 69.000000 56.0000000 153.666667  not_sel G_87
## 76 28.0000000 64.000000 68.0000000 160.000000  not_sel G_12
## 77 17.3333333 55.000000 88.0000000 160.333333  not_sel G_76
## 78 53.0000000 19.000000 89.0000000 161.000000  not_sel G_80
## 79 23.0000000 62.000000 84.0000000 169.000000  not_sel G_95
## 80 25.0000000 70.000000 75.0000000 170.000000  not_sel G_18
## 81 84.0000000 25.333333 61.0000000 170.333333  not_sel G_22
## 82 70.0000000 25.333333 77.0000000 172.333333  not_sel G_78
## 83 86.0000000 63.000000 23.6666667 172.666667  not_sel G_71
## 84 74.0000000 70.000000 30.0000000 174.000000  not_sel G_96
## 85 86.0000000 17.000000 73.0000000 176.000000  not_sel G_83
## 86 90.0000000 76.000000 10.3333333 176.333333  not_sel G_24
## 87 30.0000000 76.000000 74.0000000 180.000000  not_sel G_23
## 88 81.0000000 76.000000 24.0000000 181.000000  not_sel G_66
## 89 79.0000000 12.333333 90.0000000 181.333333  not_sel G_56
## 90 72.0000000 24.000000 86.0000000 182.000000  not_sel G_86
## 91 30.0000000 76.000000 78.0000000 184.000000  not_sel G_37
## 92 90.0000000 25.333333 69.0000000 184.333333  not_sel G_59
## 93 90.0000000 14.666667 81.0000000 185.666667  not_sel  G_7
## 94 25.3333333 76.000000 85.0000000 186.333333  not_sel G_79
## 95 30.0000000 64.000000 93.0000000 187.000000  not_sel G_15
## 96 90.0000000 76.000000 31.0000000 197.000000  not_sel G_69
## 97 30.0000000 76.000000 93.0000000 199.000000  not_sel G_62
gen.sel = df[which(df$sel=="selected"),"gen"]

BPSI plot - w 3,1,1

Highlighting the ranking of probability of superior performance per trait of the selected families.

library(knitr)

plot(BPSI_soy,category = "Ranks")
Ranking of probability of superior performance per trait - w (3,1,1)

Ranking of probability of superior performance per trait - w (3,1,1)

plot(BPSI_soy,category = "BPSI")
## Warning: Could not calculate the predicate for layer 3; ignored
Selected families based on BPSI rank - w (3,1,1)

Selected families based on BPSI rank - w (3,1,1)

library(kableExtra)
BPSI_soy$BPSI |> kable(format = "html", booktabs = TRUE, longtable=TRUE)
GY PH PL bpsi sel gen
G_20 0.3333333 4.000000 4.0000000 8.333333 selected G_20
G_38 3.0000000 2.000000 11.6666667 16.666667 selected G_38
G_3 8.0000000 6.000000 2.6666667 16.666667 selected G_3
G_4 8.6666667 3.000000 15.0000000 26.666667 selected G_4
G_45 4.6666667 5.000000 18.0000000 27.666667 selected G_45
G_13 7.0000000 8.000000 16.3333333 31.333333 selected G_13
G_16 29.0000000 12.000000 3.3333333 44.333333 selected G_16
G_1 12.6666667 18.000000 14.0000000 44.666667 selected G_1
G_33 19.0000000 14.333333 12.0000000 45.333333 selected G_33
G_14 27.0000000 4.666667 19.0000000 50.666667 selected G_14
G_32 2.0000000 49.000000 1.0000000 52.000000 not_sel G_32
G_21 15.0000000 38.000000 3.0000000 56.000000 not_sel G_21
G_28 17.0000000 11.666667 28.0000000 56.666667 not_sel G_28
G_31 6.6666667 17.000000 33.0000000 56.666667 not_sel G_31
G_54 57.0000000 1.000000 1.6666667 59.666667 not_sel G_54
G_50 1.3333333 52.000000 7.0000000 60.333333 not_sel G_50
G_93 5.0000000 47.000000 8.6666667 60.666667 not_sel G_93
G_64 48.0000000 2.333333 11.0000000 61.333333 not_sel G_64
G_60 23.0000000 25.000000 15.3333333 63.333333 not_sel G_60
G_40 13.0000000 42.000000 9.0000000 64.000000 not_sel G_40
G_68 13.6666667 39.000000 13.0000000 65.666667 not_sel G_68
G_49 34.0000000 24.000000 7.6666667 65.666667 not_sel G_49
G_39 37.0000000 12.000000 17.0000000 66.000000 not_sel G_39
G_77 24.0000000 34.000000 12.3333333 70.333333 not_sel G_77
G_8 36.0000000 29.000000 5.3333333 70.333333 not_sel G_8
G_47 21.0000000 25.333333 30.0000000 76.333333 not_sel G_47
G_36 18.0000000 16.000000 47.0000000 81.000000 not_sel G_36
G_44 9.0000000 18.666667 55.0000000 82.666667 not_sel G_44
G_53 2.0000000 61.000000 20.0000000 83.000000 not_sel G_53
G_29 15.3333333 30.000000 39.0000000 84.333333 not_sel G_29
G_52 33.0000000 13.666667 38.0000000 84.666667 not_sel G_52
G_74 43.0000000 28.000000 14.3333333 85.333333 not_sel G_74
G_58 44.0000000 13.000000 29.0000000 86.000000 not_sel G_58
G_9 9.3333333 76.000000 2.0000000 87.333333 not_sel G_9
G_25 12.0000000 24.000000 53.0000000 89.000000 not_sel G_25
G_11 32.0000000 6.333333 51.0000000 89.333333 not_sel G_11
G_42 7.3333333 58.000000 25.0000000 90.333333 not_sel G_42
G_46 11.0000000 76.000000 13.6666667 100.666667 not_sel G_46
G_10 68.0000000 23.000000 12.0000000 103.000000 not_sel G_10
G_67 73.0000000 6.666667 24.0000000 103.666667 not_sel G_67
G_30 49.0000000 10.666667 45.0000000 104.666667 not_sel G_30
G_43 10.0000000 76.000000 19.6666667 105.666667 not_sel G_43
G_75 77.0000000 8.666667 21.0000000 106.666667 not_sel G_75
G_88 65.0000000 15.000000 26.6666667 106.666667 not_sel G_88
G_92 10.0000000 53.000000 44.0000000 107.000000 not_sel G_92
G_6 22.3333333 22.000000 63.0000000 107.333333 not_sel G_6
G_89 31.0000000 16.666667 60.0000000 107.666667 not_sel G_89
G_2 55.0000000 3.000000 50.0000000 108.000000 not_sel G_2
G_17 54.0000000 3.333333 54.0000000 111.333333 not_sel G_17
G_5 81.0000000 3.666667 29.0000000 113.666667 not_sel G_5
G_70 20.3333333 72.000000 22.0000000 114.333333 not_sel G_70
G_35 45.0000000 60.000000 16.0000000 121.000000 not_sel G_35
G_34 5.3333333 76.000000 40.0000000 121.333333 not_sel G_34
G_63 35.0000000 76.000000 10.6666667 121.666667 not_sel G_63
G_26 13.3333333 76.000000 34.0000000 123.333333 not_sel G_26
G_94 47.0000000 10.333333 66.0000000 123.333333 not_sel G_94
G_91 42.0000000 25.333333 57.0000000 124.333333 not_sel G_91
G_61 59.0000000 15.000000 52.0000000 126.000000 not_sel G_61
G_27 51.0000000 76.000000 0.3333333 127.333333 not_sel G_27
G_51 50.0000000 76.000000 2.0000000 128.000000 not_sel G_51
G_41 39.0000000 25.333333 64.0000000 128.333333 not_sel G_41
G_19 86.0000000 16.000000 27.3333333 129.333333 not_sel G_19
G_84 20.6666667 40.000000 76.0000000 136.666667 not_sel G_84
G_98 21.3333333 27.000000 90.0000000 138.333333 not_sel G_98
G_65 25.6666667 21.000000 93.0000000 139.666667 not_sel G_65
G_97 80.0000000 18.000000 42.0000000 140.000000 not_sel G_97
G_81 22.0000000 59.000000 62.0000000 143.000000 not_sel G_81
G_72 63.0000000 15.333333 65.0000000 143.333333 not_sel G_72
G_85 56.0000000 64.000000 23.3333333 143.333333 not_sel G_85
G_82 58.0000000 64.000000 26.3333333 148.333333 not_sel G_82
G_73 8.3333333 76.000000 67.0000000 151.333333 not_sel G_73
G_90 60.0000000 64.000000 27.3333333 151.333333 not_sel G_90
G_57 27.0000000 33.000000 93.0000000 153.000000 not_sel G_57
G_48 23.6666667 72.000000 58.0000000 153.666667 not_sel G_48
G_87 28.6666667 69.000000 56.0000000 153.666667 not_sel G_87
G_12 28.0000000 64.000000 68.0000000 160.000000 not_sel G_12
G_76 17.3333333 55.000000 88.0000000 160.333333 not_sel G_76
G_80 53.0000000 19.000000 89.0000000 161.000000 not_sel G_80
G_95 23.0000000 62.000000 84.0000000 169.000000 not_sel G_95
G_18 25.0000000 70.000000 75.0000000 170.000000 not_sel G_18
G_22 84.0000000 25.333333 61.0000000 170.333333 not_sel G_22
G_78 70.0000000 25.333333 77.0000000 172.333333 not_sel G_78
G_71 86.0000000 63.000000 23.6666667 172.666667 not_sel G_71
G_96 74.0000000 70.000000 30.0000000 174.000000 not_sel G_96
G_83 86.0000000 17.000000 73.0000000 176.000000 not_sel G_83
G_24 90.0000000 76.000000 10.3333333 176.333333 not_sel G_24
G_23 30.0000000 76.000000 74.0000000 180.000000 not_sel G_23
G_66 81.0000000 76.000000 24.0000000 181.000000 not_sel G_66
G_56 79.0000000 12.333333 90.0000000 181.333333 not_sel G_56
G_86 72.0000000 24.000000 86.0000000 182.000000 not_sel G_86
G_37 30.0000000 76.000000 78.0000000 184.000000 not_sel G_37
G_59 90.0000000 25.333333 69.0000000 184.333333 not_sel G_59
G_7 90.0000000 14.666667 81.0000000 185.666667 not_sel G_7
G_79 25.3333333 76.000000 85.0000000 186.333333 not_sel G_79
G_15 30.0000000 64.000000 93.0000000 187.000000 not_sel G_15
G_69 90.0000000 76.000000 31.0000000 197.000000 not_sel G_69
G_62 30.0000000 76.000000 93.0000000 199.000000 not_sel G_62